clear;
clc;
tstart=1; %start time
tmax=200; %max time of simulation
time=tstart:tmax; %time vector


%nparticles = 1000;
nparticles = 1;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(1,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(1,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(1,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(1,npos);%a's for spin down at each N at current time

%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=tmax/2:100:tmax;
sigma_count = 1;
T=time;
pos_counter=0;
a=0;

avg=0;
alpha=0.005;

p_n_tmax = zeros(1,npos);

%for all particles
for p=1:nparticles
    %form the initial a vectors for time = 0. consider all the particles to be 
    %with spin up
    a = 1/2;
%     a = 1/(sqrt(2));
    a_n_u_prev(nmid_rel)=a;
    a_n_d_prev(nmid_rel)=(a)*i;
    
    
    %for all times
    for t=time
        display(t);
        pos_counter = 0;
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
        
%         alpha1=normrnd(avg,alpha);
%         alpha2=normrnd(avg,alpha);
%         alpha3=normrnd(avg,alpha);

        %for all possible positions for this particular time, fill up the a
        %arrays for the current time
%         B1=cos(alpha3) + i*sin(alpha3);
%         B2=cos(alpha1-i*alpha2) + i*sin(alpha1 - i*alpha2);
%         B3=cos(alpha1+i*alpha2) + i*sin(alpha1 + i*alpha2);
%         B4 = cos(alpha3) - i*sin(alpha3);
        
        for pos_counter=nstart_rel:nend_rel
              %Below is for the normal quantum walk of the paper without noise.  
            a_n_u_cur(pos_counter) =  (1/sqrt(2))*( a_n_d_prev(pos_counter-1) +  a_n_u_prev(pos_counter-1) );
            a_n_d_cur(pos_counter)= (1/sqrt(2))*( a_n_u_prev(pos_counter+1) - a_n_d_prev(pos_counter+1) );

%             a_n_u_cur(pos_counter) = (1/sqrt(2))* ( a_n_u_prev(pos_counter-1) * (B1 + B3) + a_n_d_prev(pos_counter-1) *(B2 + B4) );
%             a_n_d_cur(pos_counter)=  (1/sqrt(2))* (a_n_u_prev(pos_counter+1) *(B1-B3) - a_n_d_prev(pos_counter+1) *(B2-B4));

            p_n_tmax(pos_counter) = (abs(a_n_u_cur(pos_counter)))^2 + (abs(a_n_d_cur(pos_counter)))^2; 
            
        end;
            sigma(sigma_count) = sqrt( sum(pos.^2 .* p_n_tmax) - (sum(pos.*p_n_tmax)).^2 );
            sigma_count = sigma_count + 1;
            
        if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(1,npos);
            a_n_d_cur = zeros(1,npos);
        end;
        
    end;
    
end;

figure(1)
plot(pos-nmid_rel,p_n_tmax);
grid on;

figure(2)
plot(T,sigma);
grid on;

